#ifndef UHM_SOLVER_OPS_SOLVER_HXX_
#define UHM_SOLVER_OPS_SOLVER_HXX_

namespace uhm {

  template<class S_> 
  Int_ op_uhm_set_matrices(S_ *s, typename S_::_Tree_::Tnd_* me) {
    typename S_::_Tree_::_Data_ e = me->get_data();
    e->uhm().set(s->get_type(), e->factor().dof(), e->schur().dof(), 
                 s->get_rhs(), s->get_blk());
    return true;
  }

  template<class S_> 
  Int_ op_uhm_create_matrices(S_ *s, typename S_::_Tree_::Tnd_* me) {
    typename S_::_Mesh_::_Uhm_ &uhm = me->get_data()->uhm();
    uhm.create_flat_without_buffer();
    uhm.create_flat_buffer(); // do not chreate ABR                                        
    uhm.create_hier();
    return true;
  }

  template<class S_> 
  Int_ op_uhm_set_and_create_matrices(S_ *s, typename S_::_Tree_::Tnd_* me) {
    typename S_::_Tree_::_Data_ e = me->get_data();
    e->uhm().set(s->get_type(), e->factor().dof(), e->schur().dof(),
                 s->get_rhs(), s->get_blk());

    typename S_::_Mesh_::_Uhm_ &uhm = me->get_data()->uhm();
    uhm.create_flat_without_buffer();
    uhm.create_flat_buffer(); // do not chreate ABR                                        
    uhm.create_hier();
    return true;
  }

  template<class S_> 
  Int_ op_uhm_set_offset(S_ *s, typename S_::_Tree_::Tnd_* me) {
    typename S_::_Tree_::_Data_ e = me->get_data();
    s->set_offset(e->set_nodes_offset(s->get_offset()));
    return true;
  }

  template<class S_> 
  Int_ op_uhm_direct_decompose(S_ *s, typename S_::_Tree_::Tnd_* me) {
    typename S_::_Mesh_::_Uhm_ &uhm = me->get_data()->uhm();
    UHM_ERROR(op_uhm_merge_A<S_>(s, me), UHM_ERROR_SOLVER_MERGE);
    UHM_ERROR(uhm.decompose(), UHM_ERROR_DIRECT_DECOMPOSE);
    return true;
  }

  template<class S_> 
  Int_ op_uhm_up_direct_decompose(S_ *s, typename S_::_Tree_::Tnd_* me) {

    Bool_ reusable = true;

    typename S_::_Tree_::_Data_ elt = me->get_data();
    typename S_::_Mesh_::_Uhm_ &uhm = elt->uhm();

    if (!me->is_leaf()) {
      for (Int_ i=0;i<me->get_n_children();++i) {
        typename S_::_Tree_::_Data_ c_elt = me->get_child(i)->get_data();
        reusable &= (c_elt->get_marker());
        c_elt->set_marker(S_::_Mesh_::ELEMENT_MARKER_ZERO);
      }
    }

    // This entry actually shares buffer with uhm. So create now
    typename S_::_Cache_::Entry_ in(elt->factor(), elt->schur(), uhm), *out;
    Int_ cached = s->work().pull(in, &out);

    if (reusable && cached) {
      if (uhm.ABR().is_flat_buffer_null())
        uhm.ABR().create_flat_buffer();
      
      copy_A(out->uhm(), uhm);

      elt->set_marker(S_::_Mesh_::ELEMENT_REUSED);

      Double_ flop; uhm.decompose(flop);

      elt->set_flop(flop);

    } else {

      UHM_ERROR(op_uhm_merge_A<S_>(s, me), UHM_ERROR_SOLVER_MERGE);
      UHM_ERROR(uhm.decompose(), UHM_ERROR_DIRECT_DECOMPOSE);

      elt->set_flop(0.0);      

      // Caching matrices 
      if (me->get_gen() <= s->get_caching()) {
        cached = s->cache().push(in, &out);
        out->uhm().create_flat_without_buffer();
        out->uhm().create_flat_buffer(true); 
        
        copy_A(uhm, out->uhm());
      }

    }

    return true;
  }

  template<class S_> 
  Int_ op_uhm_direct_solve_0_x(S_ *s, typename S_::_Tree_::Tnd_* me) {
    typename S_::_Mesh_::_Uhm_ &uhm = me->get_data()->uhm();
    UHM_ERROR(op_uhm_merge_B<S_>(s, me), UHM_ERROR_SOLVER_MERGE);
    return true;
  }

  template<class S_> 
  Int_ op_uhm_direct_solve_1_x(S_ *s, typename S_::_Tree_::Tnd_* me) {
    typename S_::_Mesh_::_Uhm_ &uhm = me->get_data()->uhm();
    UHM_ERROR(op_uhm_merge_X<S_>(s, me), UHM_ERROR_SOLVER_MERGE);
    UHM_ERROR(uhm.solve_1_x(), UHM_ERROR_DIRECT_SOLVE_ONE);
    return true;
  }

  template<class S_> 
  Int_ op_uhm_direct_solve_2_x(S_ *s, typename S_::_Tree_::Tnd_* me) {
    typename S_::_Mesh_::_Uhm_ &uhm = me->get_data()->uhm();
    UHM_ERROR(uhm.solve_2_x(), UHM_ERROR_DIRECT_SOLVE_TWO);
    UHM_ERROR(op_uhm_branch_X<S_>(s, me), UHM_ERROR_SOLVER_BRANCH);
    return true;
  }

  template<class S_> 
  Int_ op_uhm_direct_solve_0_r(S_ *s, typename S_::_Tree_::Tnd_* me) {
    typename S_::_Mesh_::_Uhm_ &uhm = me->get_data()->uhm();
    UHM_ERROR(op_uhm_copy_R<S_>(s, me), UHM_ERROR_SOLVER_COPY);
    return true;
  }

  template<class S_> 
  Int_ op_uhm_direct_solve_1_r(S_ *s, typename S_::_Tree_::Tnd_* me) {
    typename S_::_Mesh_::_Uhm_ &uhm = me->get_data()->uhm();
    UHM_ERROR(op_uhm_merge_R<S_>(s, me), UHM_ERROR_SOLVER_MERGE);
    UHM_ERROR(uhm.solve_1_r(), UHM_ERROR_DIRECT_SOLVE_ONE);
    return true;
  }

  template<class S_> 
  Int_ op_uhm_direct_solve_2_r(S_ *s, typename S_::_Tree_::Tnd_* me) {
    typename S_::_Mesh_::_Uhm_ &uhm = me->get_data()->uhm();
    UHM_ERROR(uhm.solve_2_r(), UHM_ERROR_DIRECT_SOLVE_TWO);
    UHM_ERROR(op_uhm_branch_R<S_>(s, me), UHM_ERROR_SOLVER_BRANCH);
    return true;
  }

}


#endif
